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ABSTRACT 

We report the results of spectroscopic mapping observations carried out to- 
ward protostellar outflows in the BHR71, L1157, L1448, NGC 2071, and VIA 
1623 molecular regions using the Infrared Spectrograph (IRS) of the Spitzer Space 
Telescope. These observations, covering the 5.2 — 37 /im spectral region, provide 
detailed maps of the 8 lowest pure rotational lines of molecular hydrogen and of 
the [SI] 25.25 /xm and [Fell] 26.0 /im fine structure lines. The molecular hydro- 
gen lines, believed to account for a large fraction of the radiative cooling from 
warm molecular gas that has been heated by a non-dissociative shock, allow the 
energetics of the outflows to be elucidated. Within the regions mapped towards 
these 5 outflow sources, total H2 luminosities ranging from 0.02 to 0.75 Lq were 
inferred for the sum of the 8 lowest pure rotational transitions. By contrast, the 
much weaker [Fell] 26.0 /im fine structure transition traces faster, dissociative 
shocks; here, only a small fraction of the fast shock luminosity emerges as line 
radiation that can be detected with Spitzer /YRS>. 

Subject headings: ISM: Molecules — ISM: Abundances — ISM: Clouds - molec- 
ular processes - shock waves 
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1. Introduction 

It has long been recognized that the formation of stars in molecular clouds is often 
accompanied by supersonic outflows (e.g. Bally & Lada 1983). Where protostellar outflows 
strike nearby molecular material, shock waves can be created; these compress and heat the 
gas, and can modify its composition. Interstellar shocks have been the subject of extensive 
study, both observational and theoretical (summarized, for example, in the review of Draine 
& McKee 1993). Slow shocks - propagating at speeds less than ~ 50kms~^ (in clouds with 
typical magnetic fields) - leave the gas predominantly molecular in their wake, although the 
chemical composition can be modified by two processes: (1) the sputtering of dust grains 
and their icy mantles, which releases material into the gas phase; and (2) the enhancement of 
gas-phase chemical reactions that possess activation energy barriers; these may be negligibly 
slow at the typical temperatures of molecular clouds but rapid at the elevated temperatures 
behind shock waves. Faster shocks can lead to the coUisional dissociation of molecules and 
the ionization of the resultant atoms. 

The Infrared Spectrograph (IRS) of the Spitzer Space Telescope provides a valuable tool 
for observing emissions from shocked interstellar gas. With its coverage of 5.2—37 /zm spectral 
region, IRS has provided access to (1) the lowest 8 pure rotational transitions of molecular 
hydrogen (e.g. Neufeld et al. 2006a, hereafter N06) - the dominant constituent of molecular 
clouds; (2) pure rotational transitions of HD (e.g. Neufeld et al. 2006b), H2O, and OH (e.g. 
Melnick et al. 2008); (3) rovibrational transitions of CO2 and C2H2 (e.g. Sonnentrucker et 
al. 2006; Sonnentrucker, Gonzalez-Alfonso & Neufeld 2007); (4) fine structure emissions of 
several atoms and atomic ions (Neufeld et al. 2007, hereafter N07). Observations of shocked 
material at mid- infrared wavelengths offer significant advantages over ultraviolet, optical or 
even near-infrared studies: they are affected less by interstellar extinction, and they can 
probe cooler material than can be observed at shorter wavelengths. 

In this paper, we present the results of spectroscopic mapping observations carried out 
using the Spitzer/IRS) toward shocked material associated with five protostellar outflows: 
BHR71, L1157, L1448, NGC 2071, and VLA 1623. In §2 below, we discuss the properties 
of the individual sources and the rationale underlying our source selection. In §3, we de- 
scribe the observations, discuss our data reduction procedure, and present average spectra 
and spectral line maps for each source. In §4, we discuss the basic features of the II2 rota- 
tional emissions and atomic fine structure emissions that we detected; here we analyse the 
overall energetics for each source, describe the correlations among the various detected line 
emissions, and compare our H2 spectral line maps with continuum maps obtained previously 
using Spitzer's Infrared Array Camera (IRAC) and extracted by us from the Spitzer archive. 
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2. Source selection 

The sample selected for the present investigation is composed of five well-studied pro- 
tostellar outflows, namely L1448, BHR71, VIA 1623, L1157 and NGC 2071, that are known 
to contain warm molecular gas from ground-based observations of H2 vibrational emissions. 

Our sample is not intended to provide an unbiased survey of protostellar outflows. The 
first four sources listed above are all outflows driven by low-mass, low luminosity (L < 70 Lq) 
protostars that provide a particularly "clean" laboratory for studying shock physics. The 
fifth source, NGC 2071, contains a cluster of intermediate- mass protostars with a higher 
luminosity [L ~ 500 L©) and thus allows us to contrast the outflows associated with low- 
mass objects with those driven by intermediate mass protostars. The new observations of 
this source complement previous Spitzer observations of a single 1' x 1' region that have led 
recently to the first detection of spatially-extended water (Melnick et al. 2008). 

All the sources in our target list have been extensively studied at several wavelengths 
using a variety of different facilities; and all of them will be observed as part of the Her- 
schel Guaranteed Time Key Program Water in Starforming Regions with Herschel (WISH) . 
We have provided a description of NGC 2071 in a recent paper (Melnick et al. 2008); the 
properties of the other four sources are reviewed below. 

2.1. BHR71 

BHR 71 (also known as DC 297.2-2.8 or Sa 136) is a well-isolated Bok globule located 
at a distance of ~ 200 pc toward the constellation Musca (Sandqvist 1977; Bourke, Hyland, 
& Robinson 1995a; Bourke et al. 1995b). CO observations of this source reveal a highly 
coUimated bipolar outflow (Bourke et al. 1997) associated with IRAS 11590-6452, which 
was classified as a very young Class protostar with a luminosity of ~ 9 . Higher spatial 
resolution observations using ISOCAM, combined with ground-based near-infrared i^'-band 
(2.11 nm) imaging and ^^CO J = 2 ^ 1 observations, show that the IRAS source is likely 
composed of two embedded protostars, IRS 1 and IRS 2, with a projected separation of 
~ 17" (3400 AU; Bourke 2001). IRS 1 and IRS 2 each drive a molecular outflow, though the 
significantly larger outflow, and the one of greater interest here, is driven by IRS 1. 

Recent observations of the circumstellar envelope around IRS 1 and IRS 2 by Chen et 
al. (2008) using both the Australian Telescope Compact Array and Spitzer, reveal masses 
of 2.12 and 0.05 Mq, hydrogen (H+2H2) densities of 2.6 and 2.2x10^ cm~^, and hydrogen 
column densities of 9.4 and 2.6x10^^ cm~^, respectively. Fits to the spectral energy distri- 
bution of IRS 1 and IRS 2 based on these more recent data imply bolometric luminosities 
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for these objects of 13.5 and 0.5 Lq, respectively. 

The Spitzer/IRAC image of BHR 71 is shown in Figure 1. Larger-scale, multi-transition 
CO maps reveal that the BHR 71 outflow extends ± 4' from IRS 1. Excitation analysis implies 
30 < T <50 K and 10^ < n(H2) <3xl0^ cm'^ for the IRS 1 outflow. Though smaller, 
the outflow from IRS 2 appears to contain warmer gas, i.e., T > 300 K (Parise et al. 2006). 
The CO column densities associated with the outflow lobes from IRS 1 are estimated to be 
2-3x10^^ cm~^. Assuming n(C0)/n(H2) = 10~^, the inferred H2 column density toward 
each lobe is 2-3x10^^ cm~^. Bourke et al. (1997) find the observed morphology and velocity 
structure of the flow are well accounted for with a simple model of a biconical outflow with 
a semi-opening angle of 15°, in which the gas moves outward with a constant radial velocity 
(with respect to the cone apex) of ~ 28 km s~^ and that is viewed nearly perpendicular to 
its symmetry axis. The inclination of the outflow axis from the line of sight is found to be 
~84°. The flow masses within each lobe, taking into account the mass in the velocity range 
of the ambient cloud and optical depth effects of the flowing gas, are 0.3 and 1.0 Mq for 
the southern and northern lobes, respectively. The mechanical luminosity of the molecular 
outflow, derived from ^^CO (1-0) and ^^CO (1-0) observations, is ~ 0.5 Lq and its dynamical 
age is ~10^ yr (Bourke et al. 1997). 



2.2. L1157 

The well-collimated molecular outflow in LI 157 lies deeply embedded within an extended 
circumstellar envelope with which it is interacting (Gueth et al. 2003; Beltran et al. 2004). 
The driving source is believed to be a Class 0, low-luminosity (11 Lq) IRAS source, IRAS 
20386-1-6751 (Andre et al. 1993). The outflow is oriented approximately north-south and 
extends over about 5'. The distance to L1157 is assumed to share that of other clouds 
in Cepheus. Unfortunately, estimates for the distances of these clouds cover a fairly wide 
range: 200 to 450 pc (e.g. Kun 1998, Looney et al. 2007). Though the distance appropriate 
for LI 157 is uncertain and remains a matter of debate, we adopt a value of 440 pc (Viotti 
1969). 

Figure 2 shows the Spitzer/IKAC image of LI 157. There are notable asymmetries 
between the northern and southern lobe, but the morphology is clearly S-shaped indicating 
the presence of an underlying precessing jet. Modeling of the outflow provides evidence for 3 
or 4 separate episodes of mass ejection from the driving source over a period of ~ 15,000 years 
(Bachiller et al. 2001). Observations of a number of excited molecular transitions clearly 
indicate that shocks are present in both lobes. Near-infrared H2 v=l-0 S(l) -|- continuum 
imaging (Davis & Eisloffel 1995) reveal a number of knots of H2 emission corresponding to the 
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peaks seen in the Spitzer/\YiAC image, with the most prominent of these H2 emission peaks 
lying approximately 1.0' to the southeast and 1.2' to the northwest of IRAS 20386+6751. 
Similarly, enhanced SiO and NH3 (3,3) emission has been detected associated with the 
outflow (Mikami et al. 1992; Bachiller, Martin-Pintado, & Fuente 1993), consistent with 
elevated temperatures and an abundance increase tied to the passage of a shock. Moreover, 
the spatial coincidence of the SiO and NH3 (3,3) emission with that of the H2 v=l-0 S(l) 
emission confirm that these emission lines arise in shocked material along the outflow. 

Multi-transition CO observations of the blueshifted southern lobe suggest that the bulk 
of the outflow consists of moderately dense gas of n(H2) =(l-3)xl0^ cm~^ heated to 
^kin = 50-170 K (Hirano & Taniguchi 2001). The intense thermal emission lines of SiO 
and CS preferentially sample more compact regions within the CO-defined outflow lobe with 
densities, n(H2), of 1-5x10^ cm"^ (see Mikami et al. 1992; Nisini et al. 2007). Bachiller 
& Perez Gutierrez (1997) have obtained estimates of the column density for a number of 
molecular species toward the driving source and two positions in the southern lobe. Using 
the reported column densities for ^^CO (1-0, 2-1) and assuming an abundance relative to 
H2 of 1.1x10"^ we obtain a total H2 column density of 4.7x10^^ cm~^ toward the center of 
the outflow (a = 20^ 39"" 06^^19; 5 = +68° 02' 15^'9, J2000) and 1.3x10^2 cm'^ Aa = 20", 
A6 = —60" and 1.1 xlO^^ cm"^ Aa = 35", A6 = —95" offset from the center position. 
The outflow mass and kinetic energy (not corrected for projection effects) derived from these 
CO observations are 0.36 Mq and 1.8x10^^ erg, respectively, for the southern lobe and 0.26 
Mq and 1.9x10^^ erg, respectively, for the northern lobe. 



2.3. L1448 



In the vicinity of Lynds 1448 there are a number of molecular outflows that have 
been detected via a va riety of observational techniques in the millimeter and infrared (e.g., 
Wolf-Chase et al.ll2000[). One of the more spectacular outflows is associated with L1448-mm 
(L1448C) (IBachiller et al.lll99ll ). This flow, associated with a Class source, has been de- 
tected at a wide range of wavelengths and exhibits strong signatures of high velocity gas, 
shock heating, and also shock chemistry. The flow was first detected in CO emission with the 
high est velocities 70 k m s~^) displaying a high degree of coUimation in emission morphol- 
ogy (IBachiller et al.lll990r). The high velocity CO gas is coincident with strong H2 vibrational 
emission jBallv et al.lll993h. with H2 e mission appearing illuminated on the leading edges of 
a bipolar cavity (IBachiller et al. 1119951 ). The shocked gas is a ssociated with an active chem- 



istry with elevat ed abundances of SiO (IBachiller et al.lll99ll ). CH3OH (jjimenez-Serra et al. 
2005h . and H2O jNisini et al.lll999h . 
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In Figure 3 we provide a finder chart illustrating our Spitzer IRS map coverage of the 
L1448-mm flow and also delineate notable sources in the region. L1448-mm or L1448C is 
itself a binary system, as identified in Spitzer images, and the flow origi i iates in the northern 
source, labelled as L1448C(N) in the nomenclature of |j0rgensen et al.l (120061 ) . At the tip of 
the northern blue-shifte d lobe lies another Class protostar, L1448-IRS3 (L1448-N). This 
source is also a binary (ILooney et al.ll2000l ) with outflow emission that corresponds to the 
bright spots labelled as IRS3A and 3B on the IRAC band 2 image shown in Figure 3. 

Below we provide some specifics regarding fl ow properties. The total mechanical lum i- 



nosity of the L1448-mm flow is ~ 1.5 — 3.5 'Lrr^ flBachi 



The mass loss rate has some uncertainty. iNisini et al 



ler et al. I Il995l : iNisini et all 120001 ). 
(12000|) estimate ~ 1 x lO'^ M 



yr ^ from far-IR cooling lines. Recent Spitzer data imply 



from the atomic and molecular flows, respectively ( jPionatos et al.ll2009). Th e total H2 v 





10"^ yr"^ and 10"^ Mq yr"^ 



S(l) luminosity from the L1448-mm flow is ~ 0.19 Lq (IBally et al 



1993( 1. A comparison 



of the observed high-J CO, H2O, and H2 emission detected by ISO suggests that the flow 



is consistent with arising in low velocity C-shocks (v < 20 km s ^; INisini et al.l 120001). but 
at th e tip of the bow sho ck higher velocities (~ 100 km s~^) are required (jPionatos et al. 
2009h . INisini et al.l J200o[ ) estimate that Ay ~ 6 - 9 mag towards IRS3, Ay ~ 4 - 6 mag 
towards the tip of the L14 48-mm flow, and ~ 4.7 — 7 mag in the southern lobe. The distance 
of the source is ~ 250 pc (lEnoch et al.ll2006l ). 



2.4. VLA1623 



The VLA1623 outflow wa s discovered by CO observations of cm-wave point sources 



coincident with p Oph core A (IMontmerle and Andrei Il988l : lAndre et al. Ill990 



Loren et al. 



I990I ). The molecular data revealed a remarkably collimated outflow, which was notable 
for the clear lack of infrared point sources in the region. The detection of cold, centrally- 
concentrated s ub-mm continuum e mission led this source to be known as the prototypical 



Class source fl Andre et al. 1 1 19931) . Studies of the outflow itself find ~0.1 entrained in 



the CO outflow JPent et al.lll995[ ) with an e stimated wind rnass loss rate of M,„ ~ 10 







M 

yr"^ and mechanical luminosity of ~ 1.9 L0( lBontemps et al.lll996l : lAndre et al. Ill990l ). The 
CO flow is coincident with a series of knots seen in images of H2 vibrational emission, with a 
total luminosity of LhJI-0 S(1)] ~ 0.04 JPent et al.lll995l : ICaratti o Garatti et al.ll2006h . 



In Figure 5 we show a finder chart showing the VLA1623 outflow as delineated by thes e 



known knots of H2 emission. The "H" sources follow the nomenc lature of iPent et al.l (1995) 



but we also list the A (H5) and B (H4) sources as denoted by iPavis and Eisloffell f 



Additional sources in the field have been found by a deeper H2 survey by iGomez et al. 



1995). 
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( 120031 ) ■ We note that source H5 is also known as HH313 an d that nearby in the fi eld is 
the refiection nebula associated with the Class I source GSS30 fiGrasdalen et al.lll973l ). The 
rectangles show the area mapped in our Spitzer IRS observations, which cover the northwest 
red-shifted lobe of the fiow encompassing from VLA1623 extending to knots H4 and H5. 



Davis et al.l (119991 ) obtained high spectral resolution observations of H2 S and Q vibra- 



tional lines with a focus on knot A (HH313, H4). The resolved lines allowed an exploration 
of the ortho/para ratio of the hot shocked H2 gas and they derived an o/p ratio of 3:1, 
consisten t with thermal equi librium at an estimated gas temperature in excess of 2000 K 
(see also lEisloffel et al.ll2000l ). In addition, through the use of line ratios between lines with 
similar upper state energies, they estimate a K-band extinction of = 3±1. Recent VLBA 
parallax measureme nts of young stars associated with p Oph A find a distance of 120 pc 
dloinard etaDboOsh . 



3. Observations and results 

The five outfiow sources discussed in §2 above were observed with the Short-Low (SL), 
Short-High (SH) and Long- High (LH) modules of the IRS instrument, providing coverage of 
the entire 5.2 - 37 /im spectral region accessible to IRS at the maximum spectral resolving 
power available: A/AA ~ 60 — 127 for SL and A/AA ~ 600 for SH and LH. Regions of size 
~ 1' X 1' were mapped by stepping the IRS slit perpendicular to its long axis by one-half its 
width and, in the case of SH and LH, parallel to its long axis by 4/5 (SL) or 1/5 (LH) of its 
length. In each object, several such ~ 1' x 1' regions were arranged along the outfiow axis 
to provide optimal coverage of the interaction of the outfiow and the surrounding gas. In 
the case of the SH and LH modules, we obtained spectra of offset regions that are devoid of 
strong molecular emission; these were used for background subtraction. Table 1 summarizes 
the observational details for each object. With the exception of a single ~ 1' x 1' region 
mapped previously towards each of BHR71 and NGC 2071, the observations were carried 
out in Cycle 4 of the Spitzer mission 

The data were processed at the Spitzer Science Center (SSC), using version 17.2 or 
18.0 of the processing pipeline. They were then reduced further using the SMART software 
package (Higdon et al. 2004), supplemented by additional routines (Neufeld et al. 2006, 
2007), which provide for the removal of bad pixels in the LH and SH data, the calibration 
of fluxes obtained for extended sources, the extraction of individual spectra at each position 
sampled by the IRS slit, the creation of spectral-line maps from the set of extracted spectra, 
the removal of striping resulting from variations in response from one pixel to another, and - 
in the case of the SH and LH modules - the subtraction of a background spectrum obtained 
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at an offset position devoid of strong molecular emissions. 

In Figures 6 - 10, we present, for each object, an example spectrum obtained by averaging 
the data over a simulated Gaussian beam of radius 25". The beam center positions, specified 
in Table 1, were chosen to encompass regions of strong line emission but weak or moderate 
continuum emission. The wavelength range plotted in Figure 6 - 10 is limited to 5.2 - 31 fim; 
the quality and reliability of the spectra diminish rapidly longward of 31 fim, particularly 
in objects that show a strong continuum flux. In the SL modules, each example spectrum 
shows a clear detection of the H2 S(2) - S(7) pure rotational lines; in the SH module, H2 
S(l) and S(2) are invariably observed with large line-to-continuum and signal-to-noise ratios; 
while in the LH module, H2 S(0) is typically detected, along with the ^Pi — ^ P2 25.25 /xm 
fine structure transition of atomic sulphur and the ^.07/2 — ^ -D9/2 transition of singly-ionized 
iron. 

Spectral line maps are shown in Figures 11-15. The effective angular resolution varies 
from ~ 3" (HPBW) for the SL module to ~ 10" for the LH module (see Table 5 of N06). 
For each object, we present maps of the H2 S(0), H2 S(l), H2 S(2), H2 S(3), H2 S(5), [Fell] 
25.99 fim and [SI] 25.25 /im lines, together with the sum of the H2 S(0) - S(7) emission lines 
and the archival 3.6 fim, 4.5 fim, 5.8 fim, and 8 fim IRAC maps. The H2 S(4), S(6) and 
S(7) line maps - although somewhat noisier and not shown in Figures 11-15 - all exhibit a 
morphology nearly identical to those of the H2 S(3) and S(5) maps. In the present paper, 
we consider only those spectral lines that can be mapped reliably over most of the region 
sampled by our observations. Additional lines of H2O, OH and HD - typically much weaker 
than those considered here - have been detected in parts of some our target objects (e.g. 
Melnick et al. 2008), and will be discussed further in upcoming papers. 

Each map is normalized relative to the reference intensities listed in Table 2. For the 
spectral line maps, these reference intensities are not the maximum intensities but instead 
the values that are exceeded in only 0.1% of the mapped aresQ; they are therefore somewhat 
smaller than the maximum measured intensities. For the continuum maps, the reference 
intensities were chosen to bring out the low intensity regions within the maps; they are often 
much smaller than the maximum intensities. In the spectral line maps of BHR71, a 25" 
region near the map center has been excised from observations with the LH module; here, 
the continuum intensity is too great to permit a reliable determination of the intensities of 
lines observed with the LH module. For each object, the map center is listed in Table 1, and 
the tick marks indicate offsets relative to that position, small tick marks being separated by 
10" and large tick marks by 50". North is up and East is to the left. The green circle plotted 



This choice of color scale allows the lower intensity emission to be discerned more readily. 
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on the H2 S(5) map indicates the simulated Gaussian beam for which the example spectra 
in Figures 6-10 were obtained. 

By integrating the line intensities over the entire mapped region for each object - but 
excluding the excised region in BHR71 for those lines detected with the LH module - we 
obtained the total line luminosities listed in Table 3. These assume the distance estimates 
discussed in §2 and tabulated in Table 1, and assume the line emission to be isotropic, as 
expected for optically thin lines. 

4. Discussion 

4.1. H2 rotational emissions 

The pure rotational transitions of H2 are optically thin, and thus the level populations 
are readily calculable. Our observations of the S(0) to S(7) rotational emissions provide 
estimates of the J = 2 through J = 9 column densities. In this paper, we consider only 
the total H2 luminosities for the entire mapped region in each object. We defer a detailed 
discussion of spatial variations in the line luminosities and line ratios to a future paper, 
(Nisini et al. in preparation), in which we will also discuss existing spectral line maps of 
H2 vibrational emissions. In Figure 16, we present rotational diagrams based upon the H2 
luminosities given in Table 3, corrected for extinction using the extinction estimates given 
in Table 1. Here, we plot the logarithm of the total number of molecules in each state, 
divided by the degeneracy of that state, as a function of the energy. Diamonds indicate the 
observed values, and the solid lines are results obtained from the fitting procedure described 
below. While material in thermal equilibrium at a single temperature would exhibit a linear 
behavior in such a plot, the observed rotational diagrams - like those obtained for outflow 
sources by N06 and by Maret et al. (2009), and for supernova remnants by N07 - depart 
from that behavior in two respects. First, there is a positive curvature to the rotational 
diagram, indicative of the presence of multiple temperature components. Second, except 
toward L1448, there is a characteristic zigzag behavior, most notably for the states of lowest 
energy, in which the states of even rotational quantum number J lie systematically above 
those of odd J. As in N06, N07, and Maret al. (2009), we attribute this second behavior 
to an ortho-to-para ratio - smaller than the equilibrium value of 3 - that is the relic of an 
earlier period when the gas was considerably cooler than it is now. 

We have modeled the H2 rotational diagrams in Figure 16, using a method similar to 
that described by Neufeld and Yuan (2008; hereafter NY08). As in NY08 (their eqn. 3), 
we assume a power law distribution of gas temperatures with the number of molecules at 
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temperature between T and T + dT assumed proportional to T . We assume a temperature 
range extending from T^j^ = 100 K to T^ax = 4000 K. Solving the equations of statistical 
equilibrium for the H2 level populations, and adopting the molecular data cited in NY08, we 
obtained predictions for the pure rotational line luminosities as a function of the power law 
index, b, the total mass of hydrogen molecules warmer than Tmin, Mmi^ Tmin), and the H2 
density, n(B.2). The excitation of pure rotational emissions is dominated by collisions with 
H2 (see NY08), in contrast to the case of rovibrational emissions, where excitation by atomic 
hydrogen presents an additional complication that we will discuss in a future paper (Nisini 
et al. 2009, in preparation). 

We have extended the treatment of NY08 by considering the H2 ortho-to-para ratio 
along with b, Mh2(> ^mm) and n{B.2). Here, we introduce two additional parameters: the 
initial ortho-to-para ratio, OPRq, and the time period, r, for which the gas temperature has 
been elevated. If reactive collisions with atomic hydrogen are the dominant mechanism for 
the interconversion of ortho- and para-H2, then the ortho-to-para ratio after time r is given 
by Nofl 

QP^(^) ^ QP^O p-n(H) r [fcpo+fcop] I QP^LTE / , _ „(H) r [fepo+fcop] A (^) 

1 + OPR(r) 1 + OPRo 1 + OPRlte V J' ^ ' 

where ?t,(H) is the atomic hydrogen density, kpo is the rate coefficient for para-to-ortho 
conversion, estimated as 8 x 10^^^ exp(— 3900/T) cm'^ s~^ by Schofield et al. (1967), based 
upon the laboratory experiments of (Schulz & Le Roy 1965), and kop ~ is the rate 

coefficient for ortho-to-para conversion. Thanks to the strong dependence of kpo upon the 
gas temperature, the ortho-to-para ratio is similarly temperature-dependent, the warmer gas 
components being closer to equilbrium than the cooler ones. This effect is evident in Figure 
16, the degree of zigzag being least for the highest energy states that are excited primarily 
in the warmest gas. 

We have adjusted the five parameters b, M(H2), ■n,(H2), OPRq, and ?T.(H)r in order to 
obtain the best fit to the data. Here, we minimize the the S(l) through S(7) rotational 

lines, assuming the same fractional uncertainty for each data point. We exclude the H2 S(0) 
line fiux from the analysis, because its value is rather uncertain, particularly in regions of 
strong continuum intensity where the line-to-continuum ratio becomes small. The best fit 
parameters are given in Table 3. 



^Recognizing that gas in the 100 - 300 K can contribute significantly to the pure rotational line emission 
that we detect with Spitzer, we have reduced Tmm from 300 K in NY08 to 100 K in the present work. 

■^Note that the expression given by NOG erroneously omits fcop from the term [fcpo + ^op], resulting in a 
25% error. 
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The inferred masses of warni0 H2 (i.e. at T > 100 K) range from 1.3 x 1O~^M0 in VLA 
1623 to O.I8M0 in NGC 2071. These values exclude the mass contribution from helium and 
apply only to the regions mapped by IRS. The best-fit power-law indices, b, range from 2.2 
to 3.3. These values lie somewhat below the typical values determined by NY08 in the SNR 
IC443C, and somewhat below the value b = 3.8 expected for a paraboloidal bow shock that 
is dissociative at its apex. The best fit H2 densities range from 3000 to 7000 cm~'^. It should 
be noted, however, that the range of acceptable values for b and n(H2) can be large, and that 
these two parameters are somewhat degenerate in the following sense: increasing b tends to 
decrease the predicted relative strength for the higher J transitions, whilst decreasing n{B.2) 
has a similar effect (because the level populations become increasingly subthermal). Thus, 
an increase in one of these assumed parameters can partially compensate for a decrease in 
the other. The correlated uncertainties in b and n{B.2) are presented graphically in Figure 
17, which shows contours of in the b-n{B.2) plane. Crosses indicate the minimum (i-e. 
the best fit), while the two contours represent the 68% and 95% confidence limits on each 
parameter. With the exception of L1448, for which the data are consistent with an OPR 
of 3, the objects require an initial OPR less than 3. In two objects, L1157 and NGC 2071, 
the current OPR is clearly larger for the higher-J states than for the lower-J states. In the 
context of our model for the evolution of the OPR (eqn. 1), both cases imply a lower limit 
~ 10^ [cm~'^/n(H)] yr on the time period, r, for which the gas has been warm. 



4.2. Fine structure emissions 

The Fe+ ^Dy/a D9/2 25.988 /xm and S ^Pi P2 25.249 /im fine structure lines are 
readily detectable in the LH spectra of all five outflows sources, and are strong enough to 
map. However, these transitions make a negligible contribution to the total 5.2 - 37 fim 
line luminosity. The [SI] line luminosity ranges from ~ 1 to 2 % of the total H2 S(0) - 
S(7) luminosity. For [Fell], the corresponding range is wider ~ 0.15 to 2%. In the case 
of [Fell], these percentages are much smaller than those attained in some of the supernova 
remnants (SNR) that we observed previously (N07). In the SNR 3C391, for example, the 
[Fell] 25.988 /im line flux exceeded 10% of the total H2 S(0) - S(7) flux. Furthermore, fine 
structure emissions from more highly ionized species - Ne"*", Ne"*"^, and S"*"^ for example 
- are readily detectable in SNR but not in protostellar outflow sources. These differences 
suggest that higher velocity shocks are more prevalent in supernova remnants than around 



In this paper, we use the word "warm" to describe H2 at temperatures > 100 K. This, of course, includes 
gas at lower temperatures than that which can be probed by means of traditional, ground-based studies of 
vibrationally-excited H2. 
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outflow sources. 

A visual inspection of Figures 11 - 15 suggests that the [SI] emissions are typically more 
closely correlated with the H2 emissions than are the [Fell] emissions. This impression is 
borne out by a correlation analysis. Comparing the H2 S(5) and [SI] intensities, spatial pixel 
by spatial pixel, we obtained linear correlation coefficients of 0.81, 0.78, 0.42, 0.80, and 0.85 
respectively for BHR71, L1157, L1448, NGC 2071 and VIA 1623. Comparing the H2 S(5) 
and [Fell] intensities, the corresponding results are invariably smaller: 0.55, 0.48, 0.23, 0.72, 
and 0.78. This behavior is consistent with our previous study of supernova remnants, in 
which a principal component analysis suggested that the [SI] emitting region was largely 
cospatial with the H2 emitting region, while emissions from ionized species such as Fe"*" 
originated in a separate gas component. This segregation is presumably a consequence of 
the different shock velocities needed to generate [SI], H2, and [Fell] emissions. The first 
two arise in slow non-dissociative shocks that create very little ionization, while the third 
arises in faster shocks that are dissociative and ionizing. According to theoretical models 
(e.g. Allen et al. 2008), the emission from dissociative shocks is dominated by optical and 
ultraviolet line emissions, with infrared fine structure emissions accounting for only a small 
fraction of the overall shock luminosity. 



4.3. Source energetics and comparison with IRAC maps 

The wavelength coverage of Spitzer is well suited to observing the radiative cooling of 
non-dissociative molecular shocks. The study of Kaufman and Neufeld (1996) indicates that 
over a wide range of shock parameters, pure rotational emissions from H2 are key coolants 
of non- dissociative shocks. For shocks of velocity in the range 20 to 40 km s~^ and with a 
preshock H2 density of 10^ cm~^, for example, between 45 and 69% of the shock luminosity is 
predicted to emerge in pure rotational emission from H2. Moreover, almost all of the H2 pure 
rotational emission emerges in the S(0) - S(7) lines accessible to Spitzer. The total H2 S(0) - 
S(7) luminosities listed in Table 1 range from 0.02 to 0.75 Lq. Taking these luminosities as 
lower limits for the total mechanical luminosity dissipated in the shocks, adopting 20 km s~^ 
as a typical shock velocity, Vg, and setting the mechanical luminosity equal to ^Mv"^, we 
obtain lower limits on the mass fiow rates through the shock, M, ranging from ~ 10"^ to 
~ 2 X lO-^Moyr-^ 

In some objects, H2 line emissions can contribute significantly to the intensities measured 
using broadband mid-IR photometry (e.g. Reach et al. 2006, NY08). For the supernova 
remnant IC443, NY08 found that the mid-IR spectral region covered by IRAC was dominated 
by H2 line emissions. In particular, the intensities measured in the 5.8 /im and 8 /im IRAC 
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channels (Bands 3 and 4), were accounted for almost entirely by the IRS-measured H2 S(4) 
- S(7) line intensities. Given an excitation model for H2, the 3.6 fim and 4.5 fim IRAC 
intensities could similarly be explained as resulting from H2 emissions, primarily the S(9) 
pure rotational emission for the 4.5 fim band and the v = 1 - Q(5) rovibrational line for 
the 3.6 yum band. 

For the protostellar outflow sources studied here, we can evaluate the contribution of 
the H2 S(4) - S(7) emissions to the broadband intensities measured with IRAC. In Table 3, 
we list the fraction of the intensities measured in the 5.8 fim and 8 fim IRAC bands that 
is attributable to H2 S(4) and S(5) (in the case of the 8 /xm band) or H2 S(6) and S(7) (in 
the case of the 5.8 /^m band). The contribution is clearly largest in BHR71 and L1157, for 
which the computed contribution to the 8 /im band actually exceeds 100%, a discrepancy 
that presumably reflects flux calibration errors that are within expectations. Even in these 
objects, however, the 5.8 fim band intensity is 2-3 times stronger than what can be attributed 
to H2, and in the three other objects, the contribution of H2 to either the 5.8 /im or 8.0 fim 
band is small, at least averaged over the entire mapped region. Nevertheless, the IRAC maps 
shown in Figures 11-15 can exhibit clear morphological similarities to the H2 line maps, 
even in objects such as L1448 where the overall contribution of H2 is small. This reflects the 
fact that regions of diffuse emission can be H2-dominated even if compact continuum sources 
account for most of the IRAC intensity when averaged over the entire region that we have 
mapped. Finally, we note that in the case of VLA 1623 and NGC 2071, the contributions 
listed in Table 3 reflect an important selection effect: we speciflcally avoided observing H2 
emissions in those regions known to exhibit the strongest continuum emission, because they 
would have saturated the IRS detectors. Thus, had we been able to map the entire object 
with IRS, the fractional contributions of H2 listed in Table 3 would have been even smaller, 
probably by a large factor. In summary, a comparison of the IRS and IRAC data indicates 
that the 8/im IRAC channel can be dominated by H2 emissions for a protostellar outflow, 
as it was for the supernova remnant IC443 discussed by NY08, but that the H2 contribution 
can also be negligible, as in NGC 2071 and VLA 1623. 



5. Summary and future outlook 

1) We have carried out spectroscopic mapping observations toward protostellar outflows 
in the BHR71, LI 157, L1448, NGC 2071, and VLA 1623 molecular regions using the Infrared 
Spectrograph (IRS) of the Spitzer Space Telescope. 

2) Our observations provide detailed maps of the 8 lowest pure rotational lines of molec- 
ular hydrogen and of the [SI] 25.25 /im and [Fell] 26.0 /im flne structure lines. 
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3) Within the regions mapped towards these 5 outflow sources, total H2 luminosities 
ranging from 0.02 to 0.75 Lq were inferred for the sum of the 8 lowest pure rotational 
transitions. The overall H2 luminosities were fit with a model in which an admixture of 
gas temperatures is present within the beam. For an assumed power-law distribution of gas 
temperatures, with the mass of material at temperature T to T + dT assumed proportional to 
T~^dT, the best-fit power-law indices lie in the range 2.3 to 3.3. With the exception of L1448, 
all objects exhibit a non-equilibrium H2 ortho-to-para ratio < 3. L1157 and NGC 2071 show 
evidence for an H2 ortho-to-para ratio that is an increasing function of gas temperature, as 
expected if para-to-ortho conversion possesses an activation energy barrier. 

4) As in supernova remnants studied by Neufeld et al. 2007, the [SI] 25.3 /im fine struc- 
ture emission is more strongly correlated with the H2 emissions than is the [Fell] 26.0 /im fine 
structure emission. This behavior is consistent with the suggestion that the [Fell] 26.0 /xm 
fine structure transition traces faster, dissociative shocks. 

5) In BHR71 and LI 157, H2 S(4) and S(5) line emissions account for most of the intensity 
observed in IRAC S/xm maps of the region we mapped with IRS. 

The primary emphasis in the present paper has been on object-averaged line luminosi- 
ties. However, the data presented here will allow the construction of line ratio maps that 
can probe the variation of the physical conditions within each individual object; these will 
be discussed in a series of future papers. 
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Table 1. Observational details and source properties 



BHR71 L1157 L1448 NGC 2071 VIA 1623 



Adopted distance (pc) 200 440 250 390 120 

Map center R.A. (J2000) 12h 01m 36.07s 201i 39m 5.11s 3h 25m 38.45s 5h 47m 4.05s 16h 26m 22.96s 

Map center Dec. (J2000) -65d08'50.5" 68d 02'43.9" 30d44'13.8" Od 21' 56.3" -24d23'56.2" 

Reference position R.A. " 12h 03m 40.00s 20h 38m 03.00s 31i 24m 40.00s 5h 46m 50.00s 16h 26m 11.00s 

Reference position Dec. -64d58'20.0" 68d04'0.0" 30d46'0.0" Od 11' 55.0" -24dl7'50.0" 

Start date for IRS observations 2005 Jul 02 2007 Oct 31 2008 Feb 19 2004 Oct 21 2008 Mar 27 

End date for IRS observations 2008 Feb 29 2007 Nov 01 2008 Feb 29 2008 Mar 27 2008 Mar 28 

IRS observing time (hr) 18.8 22.8 22.8 23.9 18.3 

Beam offset'' for plotted spectrum (-2, 78) (+28, -85) (34, -82) (63, 59) (-53, 51) 

Extinction estimate {Ay in mag) 2'^ 2^ 6'^ 13*^ 15"^'-^ 

'^Used for background subtraction for the SH and LH modules 
^(Aacos5, A6) relative to the map center in units of arcsec 
'^Caratti o Garatti et al. 2006 
•^Nisini et al. 2000 
^Melnick et al. 2008 
^Davis et al. 1999 
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Table 2. Reference intensities for each map*^ 





BHR71 


L1157 


L1448 


NGC 2071 


VIA 1623 


H2 S(5) 6.910 /xm 


372.1 


262.0 


256.8 


1290.0 


409.3 


H2 S(3) 9.665 /im 


239.6 


256.3 


225.7 


628.3 


183.5 


H2 S(2) 12.279 /im 


53.4 


96.7 


60.8 


300.3 


133.8 


H2 S(l) 17.035 /im 


34.5 


55.0 


49.4 


164.2 


34.1 


H2 S(0) 28.219 //m 


4.1 


8.0 


6.0 


31.8 


15.7 


H2 S(0) - S(7) total 


1225.0 


984.0 


921.8 


3912.0 


1377.0 


Fe+ 6D7/2-^D9/2 25.988 /xm 


4.5 


6.1 


20.5 


18.7 


27.0 


S ^Pi -3 P2 25.249 /im 


25.0 


12.3 


17.0 


59.5 


27.3 


IRAC Band 1 (3.6 /im cont.) 


10 


3 


5 


40 


100 


IRAC Band 2 (4.5 /xm cont.) 


20 


3 


10 


100 


100 


IRAC Band 3 (5.8 /im cont.) 


40 


10 


10 


200 


300 


IRAC Band 4 (8.0 /im cont.) 


20 


5 


10 


200 


500 



^in units of 10 ^ erg cm ^s ^ sr ^ for lines and MJy sr ^ for continuum. The errors 
in the line intensities are likely dominated by systematic effects and are expected to 
be < 25% (N06). 
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Table 3. Measured line luminosities and derived parameters 





BHR71 


L1157 


LI 448 


NGC 2071 


VLA 1623 














112 ^ \ f J O .O A. A. jJjiLL 


7 1 


10.0 






9 1 

Zi . 1. 


H2 S(6) 6.100 /im 


3.9 


8.5 


3.4 


54.4 


2.1 


H2 S(5) 6.910 /im 


11.4 


27.9 


9.2 


206.4 


5.6 


H2 S(4) 8.025 /im 


6.8 


25.9 


7.9 


98.7 


2.5 


H2 S(3) 9.665 /xm 


10.1 


42.1 


10.8 


145.7 


2.7 


H2 S(2)* 12.279 /im 


2.2 


16.2 


2.9 


79.0 


1.6 


H2 S(l) 17.035/im 


1.7 


13.0 


3.7 


58.2 


0.6 


H2 S(0) 28.219/im 


0.3 


2.3 


0.6 


12.0 


0.4 


H2 S(0) - S(7) total 


43.7 


151.2 


43.6 


748.8 


17.7 


Fe+ 6D7/2-6D9/2 25.988 /im 


0.5 


0.2 


0.7 


2.9 


0.3 


S _3p2 25.249 /im 


0.7 


1.4 


0.8 


8.8 


0.3 



Contribution of H2 to IRAC bands (percentage) 
Band 3 (5.8 /im) 33 
Band 4 (8.0 /im) 120^ 



39 
106"= 



15 
18 



9.1 
5.7 



6.5 
2.6 



H2 excitation: best-fit parameters'^ 
H2 mass above 100 K (IO'^Mq) 
Power law index, b = -dlnM/dlnT 
log(H2 density/cm~^) 
Initial H2 OPR (OPRq) 



2.5 

2.5 ±0.5 
3.8 ±0.4 
< 2.4 



49 

3.3 ±0.5 
3.8 ±0.5 

< 1.8 



2.9 ±0.5 
3.8 ±0.4 

note e 



177 
3.2 ±0.6 
3.8 ±0.4 
< 1.8 



1.3 

9 0+1.3 

^•"-•-0.8 
o c:+0.7 
■-••^-0.1 

< 1.9 



*in units of 10 ^ Lq for assumed distances of 200, 440, 250, 390, and 120 pc respectively for 
BHR71, LI 157, L1448, NGC 2071, and VLA 1623. These values do not include any extinction 
correction. The errors in the line intensities are likely dominated by systematic effects and are 
expected to be < 25% (N06). 

^The luminosity for H2 S(2) is determined from the SL module 

*^Here, the computed contribution to the 8 /im band actually exceeds 100%, a discrepancy that 
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presumably reflects flux calibration errors that are within expectations, 
and < symbols indicate 95% confldence intervals 
"^No useful constraint (consistent with 3 or any smaller value) 
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12*^02'^ 12^01"^ 30^ 

Right Ascension 



Fig. 1 - IRAC four-band grayscale image of the BHR71 region on which are superposed 
contours of integrated CO (3-2) emission (Parise et al. 2006). The positions of the two Class 
sources associated with BHR 71 are shown; IRS 1 is believed responsible for driving the 
north-south outflow, while IRS 2 is driving a much smaller east-west outflow. The rectangles 
denote the fields surveyed by Spitzer/IRS (for the LH module). 
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20^ 39"^ 20^ 20^ 39"^ 20^ 38"^ W 
Right Ascension 



Fig. 2 - IRAC greyscale image of the LI 157 region. The color stretch is shghtly exaggerated 
to emphasize channel 4 (8/im), where the extinction is the largest. Superposed on the 
Spitzer image is the CO (2-1) emission from Bachiller et al. (2001). The position of the 
driving source, IRAS 20386+6751, is shown along with the fields surveyed by Spitzer/IKS 
(LH module). 
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30°46'00" 



30°44'00" 



30°42'00" 



L1448-IRS3A 



3'^25'"50' 



L1448-mm 




1448-IRS3B 



M 



30^ 



40 

a (J2000) 

Fig. 3 - Image of Spitzer /IKAC band 2 in the region surrounding the Class source L1448. 
Prominent sources are labelled in the image. The rectangles denote the fields surveyed by 
Spitzer/IRS (LH module). 
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5h 47m 1 QS 5h 47m QQS 

Right Ascension 

Fig. 4 - Map of the outflow from NGC 2071 as traced by the H2 1-0 S(l) 2.12 fim emission 
(Eisloffel 2000). The rectangles denote the fields surveyed by Spitzer/IRS (LH module). 
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« 



. I I I I I I I t I I i_J 

16^26"'30^ 20' 
a (J2000) 

Fig. 5 - Image of Spitzer/IKAC band 2 in the region surrounding the Class source VLA1623. 
Sources labeled with an H denote knots of vibrationally excited H2. Source H5 is also known 
as HH313. The grey box shows the approximate coverage of our Spitzer/IRS maps. 
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1 2 



1 4 



1 6 



20 



22 



24 



2( 



2^ 



1 4 




20 




30 



Fig. 6 - Example spectrum for BHR71, obtained by averaging the data over a simulated 
Gaussian beam of radius 25". The beam center position, a = 12h 01m 35.75s, 6 = -65d 
07' 33.5" (J2000), was chosen to encompass a region of strong line emission but weak or 
moderate continuum emission. The wavelengths are given in units of fim, and the intensities 
in units of MJy sr~^. 




Fig. 7 - same as Fig. 6, but for LI 157. The beam center position is a = 20h 39m 10.09s, 
5 = 68d 01' 18.9" (J2000) 
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LI 44; 




5.0 



5.5 



6.0 



6.5 



'.0 



1 



20 







1 2 



1 2 



1 4 



1 6 



■.5 




22 



24 



2( 



2^ 



1 4 




20 




30 



Fig. 8 - same as Fig. 1, but for L1448. The beam center position is a = 3h 25m 41.09s, 6 
30d 42' 51.8" (J2000). 
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GC 2071 








1 2 



1 4 



200 
1 50 
1 00 
50 




300 
250 

200 

1 50 

1 00 

50 



-50 



" CN 


5(2) ^ 


=(1) ' \ 
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1 6 



1 



20 



[SI] [Fell] H2 S(0) 




20 



22 



24 



2( 
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Fig. 9 - same as Fig. 6, but for NGC 2071. The beam center position is a = 5h 47m 8.25s, 
6 = 0d 22' 55.3" (J2000). The weak features at 28.50 fim and 29.84 /im are the HD R(3) and 
H2O 725 ^ 616 transitions (discussed previously by Melnick et al. 2008) 




Fig. 10 - same as Fig. 6, but for VLA1623. The beam center position is a = 16h 26m 19.08s, 
6 = -24d 24' 47.2" (J2000) . Order mismatches are evident in the LH spectrum (bottom 
panel) . 



- 32 - 




Relative intensity 




0.0 0.2 0.4 0.6 0.8 1.0 

Fig. 11 - Spectral line maps obtained toward BHR71 for the H2 S(0), H2 S(l), H2 S(2), H2 
S(3), H2 S(5), [Fell] 25.99 /im and [SI] 25.25 fim lines, together with the sum of the H2 S(0) 
- S(7) emission lines and the archival 3.6 /im, 4.5 /im, 5.8 fim, and 8 /im IRAC maps. Each 
map is normalized relative to the reference intensities listed in Table 2. The map center 
position is a = 12h 01m 36.07s, 6 = -65d 08' 50.5" (J2000), and each small tick interval is 
10 arcsec. The green circle indicates the size and position of the 25 arcsec diameter region 
for which the example spectrum was obtained. The offset of this position relative to the line 
center is given in Table 1. 
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Fig. 12 - same as Fig. 11, but for L1157. The map center position is a = 20h 39m 5.11s, 
5 = 68d 02' 43.9" (J2000). 
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Fig. 13 - same as Fig. 11, but for L1448. The map center position is a = 3h 25m 38.45s, 
5 = 30d 44' 13.8" (J2000). 



- 35 - 



H2 s(o) 


^ H2S(i); 

%\ 


— S(2) 


— s(3) ; 


; H2 s(5) 




1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

^ 




. 3.6yu,m 

• 


. 4.5yU.m 

• 

• 


5.8yu.m 

* 

• • 

• 

• 

• 


• 

V 

• 



Relative intensity 




0.0 0.2 0.4 0.6 0.8 1.0 



Fig. 14 - same as Fig. 11, but for NGC 2071. The map center position is a = 5h 47m 4.05s, 
6 = 0d 21' 56.3" (J2000). 
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Fig. 15 - same as Fig. V. 
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Fig. 16 - Rotational diagrams for each object. The horizontal axis is the energy of each state, 
divided by Boltzmann's constant, in units of Kelvin. The vertical axis is the logarithm of 
the number of molecules in each state, divided by the product of the degeneracies associated 
with nuclear spin and with rotation. Diamonds indicate the observed values, and the solid 
lines are results obtained from the fitting procedure described in §4.1 
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Fig. 17 - Contours of in the b-n(B.2) plane (see §4.1). Crosses indicate the minimum 
(i.e. the best fit), while the two contours represent the 68% and 95% confidence limits on 
each parameter. 



